theme(
axis.title= element_text(size= 72, face= "bold"),
axis.text= element_text(size= 72, face= "bold"),
axis.text.x= element_text(
size= 72,
face= "bold",
angle= 45,
hjust= 1,
vjust= 1))
}
#Figure 1, save with width 2400 and maintain aspect ratio
df_ret= make_binned_df(winsor_beta_ratio_dji_ret, "dji_ret")
df_novol= make_binned_df(winsor_beta_ratio_dji_novol, "dji_novol")
#Add 0.5 so the lines appear correct on the plot (columns are centered at .5)
vlines_ret= c(abs(out_stats_beta_scaled[1,3]/out_stats_beta_scaled[1,2])+0.5,
abs(out_stats_beta_scaled[1,7]/out_stats_beta_scaled[1,6])+0.5,
abs(out_stats_beta_scaled[1,11]/out_stats_beta_scaled[1,10])+0.5)
vlines_novol= c(abs(out_stats_beta_scaled[2,3]/out_stats_beta_scaled[2,2])+0.5,
abs(out_stats_beta_scaled[2,7]/out_stats_beta_scaled[2,6])+0.5,
abs(out_stats_beta_scaled[2,11]/out_stats_beta_scaled[2,10])+0.5)
plot_beta_hist(df_ret, vlines_ret)
plot_beta_hist(df_novol, vlines_novol)
#Table 7
summary_table_sims_scaled= summary_table_sims
summary_table_sims_scaled[1, c(2, 4)]= summary_table_sims_scaled[1, c(2, 4)] * 10000
summary_table_sims_scaled[1, c(3, 5)]= summary_table_sims_scaled[1, c(3, 5)] * 100000000
summary_table_sims_scaled[2, c(2, 4)]= summary_table_sims_scaled[2, c(2, 4)] * 100
summary_table_sims_scaled[2, c(3, 5)]= summary_table_sims_scaled[2, c(3, 5)] * 10000
stargazer(summary_table_sims_scaled, summary=F, digits= 2)
avail_series= c("dji_ret", "dji_novol")
measures= c("measure_pos_unit", "measure_neg_unit", "measure_pess_unit")
sim_lists= list(dji_ret= sim_results_dji_ret,
dji_novol= sim_results_dji_novol)
dm_results1= data.frame(measure= character(), ret_series= character(),
dm_mae_exp1= numeric(), dm_mse_exp1= numeric(),
dm_mae_rec1= numeric(), dm_mse_rec1= numeric(),
stringsAsFactors= FALSE)
dm_results5= data.frame(measure= character(), ret_series= character(),
dm_mae_exp5= numeric(), dm_mse_exp5= numeric(),
dm_mae_rec5= numeric(), dm_mse_rec5= numeric(),
stringsAsFactors= FALSE)
oos_r2_results = data.frame(measure= character(), ret_series= character(),
oos_r2_exp_mean= numeric(), oos_r2_rec_mean= numeric(),
oos_r2_exp_sd= numeric(), oos_r2_rec_sd= numeric(),
oos_r2_exp_min= numeric(), oos_r2_exp_max= numeric(),
oos_r2_rec_min= numeric(), oos_r2_rec_max= numeric(),
stringsAsFactors = FALSE)
calc_r2_oos = function(real, bench){
real = real[!is.na(real)]
bench = bench[!is.na(bench)]
1 - mean(real^2) / mean(bench^2)
}
for (i in seq_along(avail_series)) {
current_series= avail_series[i]
sim_results= sim_lists[[current_series]]
for (j in seq_along(measures)) {
current_measure= measures[j]
model_exp= sentiment_reg(current_series, current_measure, "full",
"exp", "usrecd", sentiment_data)
model_rec= sentiment_reg(current_series, current_measure, "full",
"rec", "usrecd", sentiment_data)
actual_resid_exp= resid(model_exp)
actual_resid_rec= resid(model_rec)
nsim= length(sim_results)
oos_exp_vec= numeric(nsim)
oos_rec_vec= numeric(nsim)
for (k in seq_len(nsim)) {
sim_resid_exp= sim_results[[k]]$expansion$residuals
sim_resid_rec= sim_results[[k]]$recession$residuals
oos_exp_vec[k]= calc_r2_oos(actual_resid_exp, sim_resid_exp)
oos_rec_vec[k]= calc_r2_oos(actual_resid_rec, sim_resid_rec)
}
clean_measure= gsub("measure_|_unit", "", current_measure)
oos_r2_results= rbind(oos_r2_results,
data.frame(measure= clean_measure,
ret_series= current_series,
oos_r2_exp_mean= mean(oos_exp_vec, na.rm= TRUE),
oos_r2_rec_mean = mean(oos_rec_vec, na.rm= TRUE),
oos_r2_exp_sd = sd(oos_exp_vec, na.rm= TRUE),
oos_r2_rec_sd = sd(oos_rec_vec, na.rm= TRUE),
oos_r2_exp_min = min(oos_exp_vec, na.rm= TRUE),
oos_r2_exp_max = max(oos_exp_vec, na.rm= TRUE),
oos_r2_rec_min = min(oos_rec_vec, na.rm= TRUE),
oos_r2_rec_max = max(oos_rec_vec, na.rm= TRUE),
stringsAsFactors = FALSE))
}
}
oos_tab= oos_r2_results
oos_tab[ , sapply(oos_tab, is.numeric)]= oos_tab[ , sapply(oos_tab, is.numeric)] * 100
stargazer(oos_tab, summary=F, digits= 2, rownames=F)
stargazer(oos_tab, summary=F, digits= 3, rownames=F)
library(readr)
library(dplyr)
library(tidyr)
library(datawizard)
library(lubridate)
library(lmtest)
library(sandwich)
library(kSamples)
library(tseries)
library(forecast)
library(urca)
library(rugarch)
library(janitor)
library(stargazer)
library(ggplot2)
library(scales)
library(zoo)
#Set seed
set.seed(12191992)
#Read data
sentiment_data= read_csv("Data/sentiment_data.csv",
col_types= cols(date= col_date(format= "%Y-%m-%d")))
# sentiment_reg()
# Helper function to run sentiment predictability regressions
# Inputs:
# ret_series: character string, name of the dependent return series column.
# measure: character string, name of the sentiment measure column.
# controls: character string, determines regression specification.
#   Options: "full", matches  Garcia (2013), includes 5 lags of measure, 5
#            lags of lhs, 5 lags of lhs squared, weekday dummies
#            "reduced" weekday dummies, only 1 lag of rets and square rets
#            "only" 1 lag of measure only
# data: input data file
#
# recession: character string
#   Options: "rec", "exp"
#
# rec_ind: character string, name of recession indicator column, "usrecd"
#
# Output: lm() regression object.
#
# Example:
# model1= sentiment_reg(ret_series= "dji_ret", measure= "measure_pos_unit",
# controls= "full", recession= "rec", rec_ind= "usrecd", data= sentiment_data)
sentiment_reg= function(ret_series, measure, controls = c("full", "reduced", "only"),
recession = c("rec", "exp"), rec_ind, data) {
sentiment_data= data
controls= match.arg(controls)
recession= match.arg(recession)
if (recession == "rec") {
reg_data= sentiment_data %>%
filter(.data[[rec_ind]] == 1)
} else {
reg_data= sentiment_data %>%
filter(.data[[rec_ind]] == 0)
}
if (controls == "full") {
reg_data$weekday_factor= as.factor(reg_data$weekday)
measure_lags= paste0("lag", 1:4, "_", measure)
ret_lags= paste0("lag", 1:5, "_", ret_series)
ret_sq_lags= paste0("lag", 1:5, "_", ret_series, "_sq")
rhs= c(measure, "weekday_factor", measure_lags, ret_lags, ret_sq_lags)
} else if (controls == "reduced"){
reg_data$weekday_factor= as.factor(reg_data$weekday)
ret_lags= paste0("lag", 1:1, "_", ret_series)
ret_sq_lags= paste0("lag", 1:1, "_", ret_series, "_sq")
rhs= c(measure, "weekday_factor", ret_lags, ret_sq_lags)
}
else if (controls == "only") {
rhs= c(measure)}
regression_formula= as.formula(paste(ret_series, "~", paste(rhs, collapse= " + ")))
regression_model= lm(regression_formula, data= reg_data)
return(regression_model)
}
avail_series= c("dji_ret", "dji_novol")
measures= c("measure_pos_unit", "measure_neg_unit", "measure_pess_unit")
out_stats_beta= data.frame(ret_series= avail_series)
out_stats_mae= data.frame(series= character(),
mae_exp= numeric(), mae_rec= numeric(),
p_ttest= numeric(), p_ks= numeric(),
p_ad= numeric(), stringsAsFactors= FALSE)
out_stats_mse= data.frame(series= character(),
mse_exp= numeric(), mse_rec= numeric(),
p_ttest= numeric(), p_ks= numeric(),
p_ad= numeric(), stringsAsFactors= FALSE)
out_stats_directional= data.frame(measure= character(), ret_series= character(),
hit_exp= numeric(), p_pt_exp= numeric(),
p_ag_exp= numeric(), hit_rec= numeric(),
p_pt_rec= numeric(), p_ag_rec= numeric(),
p_ttest_hit= numeric(),
stringsAsFactors = FALSE)
row_counter= 1
for (i in 1:length(avail_series)){
for (j in 1:length(measures)){
model_exp= sentiment_reg(avail_series[i], measures[j], "full", "exp", "usrecd", sentiment_data)
model_rec= sentiment_reg(avail_series[i], measures[j], "full", "rec", "usrecd", sentiment_data)
beta_exp= coef(model_exp)[measures[j]]
beta_rec= coef(model_rec)[measures[j]]
white_exp= coeftest(model_exp, vcov. = vcovHC(model_exp, type= "HC1"))
white_rec= coeftest(model_rec, vcov. = vcovHC(model_rec, type= "HC1"))
t_exp= white_exp[measures[j], "t value"]
t_rec= white_rec[measures[j], "t value"]
resid_exp= resid(model_exp)
resid_rec= resid(model_rec)
abs_exp= abs(resid_exp)
abs_rec= abs(resid_rec)
sq_exp= resid_exp^2
sq_rec= resid_rec^2
mae_exp= mean(abs_exp, na.rm= TRUE)
mae_rec= mean(abs_rec, na.rm= TRUE)
mse_exp= mean(sq_exp, na.rm= TRUE)
mse_rec= mean(sq_rec, na.rm= TRUE)
actual_exp= model.response(model.frame(model_exp))
actual_rec= model.response(model.frame(model_rec))
forecast_exp= fitted(model_exp)
forecast_rec= fitted(model_rec)
hit_vec_exp= ifelse(sign(forecast_exp) == sign(actual_exp), 1, 0)
hit_vec_rec= ifelse(sign(forecast_rec) == sign(actual_rec), 1, 0)
mae_ttest= t.test(abs_exp, abs_rec)$p.value
mae_ks= ks.test(abs_exp, abs_rec)$p.value
mae_ad= ad.test(abs_exp, abs_rec, method= "asymptotic")$ad[1, 3]
mse_ttest= t.test(sq_exp, sq_rec)$p.value
mse_ks= ks.test(sq_exp, sq_rec)$p.value
mse_ad= ad.test(sq_exp, sq_rec, method= "asymptotic")$ad[1, 3]
pt_exp= DACTest(forecast_exp, actual_exp, test= "PT", conf.level = 0.99)
pt_rec= DACTest(forecast_rec, actual_rec, test= "PT", conf.level = 0.99)
ag_exp= DACTest(forecast_exp, actual_exp, test= "AG", conf.level = 0.99)
ag_rec= DACTest(forecast_rec, actual_rec, test= "AG", conf.level = 0.99)
p_ttest_hit= t.test(hit_vec_exp, hit_vec_rec)$p.value
measure_name= gsub("measure_|_unit", "", measures[j])
series_name = paste0(measure_name, "_", avail_series[i])
col_exp= paste0("beta_", measure_name, "_exp")
col_rec= paste0("beta_", measure_name, "_rec")
col_t_exp= paste0("t_", measure_name, "_exp")
col_t_rec= paste0("t_", measure_name, "_rec")
out_stats_beta[i, col_exp]= beta_exp
out_stats_beta[i, col_rec]= beta_rec
out_stats_beta[i, col_t_exp]= t_exp
out_stats_beta[i, col_t_rec]= t_rec
out_stats_mae[row_counter, ]= list(series_name, mae_exp, mae_rec,
mae_ttest, mae_ks, mae_ad)
out_stats_mse[row_counter, ]= list(series_name, mse_exp, mse_rec,
mse_ttest, mse_ks, mse_ad)
out_stats_directional[row_counter, ] = list(measure_name, avail_series[i],
pt_exp$DirAcc, pt_exp$p.value,
ag_exp$p.value, pt_rec$DirAcc,
pt_rec$p.value, ag_rec$p.value,
p_ttest_hit)
row_counter= row_counter + 1
}
}
#Table 1
cols_to_scale= !grepl("^t_", names(out_stats_beta)) & seq_along(out_stats_beta) != 1
out_stats_beta_scaled= out_stats_beta
out_stats_beta_scaled[1, cols_to_scale]= out_stats_beta_scaled[1, cols_to_scale] * 10000
out_stats_beta_scaled[2, cols_to_scale]= out_stats_beta_scaled[2, cols_to_scale] * 100
stargazer(out_stats_beta_scaled, summary=F, digits= 2, rownames=F)
#Table 5
out_stats_mae_scaled= out_stats_mae
out_stats_mae_scaled[1:3, 2:3]= out_stats_mae_scaled[1:3, 2:3] * 10000
out_stats_mae_scaled[4:6, 2:3]= out_stats_mae_scaled[4:6, 2:3] * 100
out_stats_mse_scaled= out_stats_mse
out_stats_mse_scaled[1:3, 2:3]= out_stats_mse_scaled[1:3, 2:3] * 100000000
out_stats_mse_scaled[4:6, 2:3]= out_stats_mse_scaled[4:6, 2:3] * 10000
stargazer(out_stats_mae_scaled, summary=F, digits= 2, rownames=F)
stargazer(out_stats_mse_scaled, summary=F, digits= 2, rownames=F)
#Table 6
stargazer(out_stats_directional, summary=F, digits= 4, rownames=F)
measures= c("measure_pos_unit", "measure_neg_unit", "measure_pess_unit")
spec= c("full", "reduced")
out_stats_r2= data.frame(measure= gsub("measure_|_unit", "", measures),
adjr2_exp_full= NA, adjr2_rec_full= NA,
r2_exp_full= NA, r2_rec_full= NA,
adjr2_exp_reduced= NA,adjr2_rec_reduced= NA,
r2_exp_reduced= NA, r2_rec_reduced= NA,
stringsAsFactors= FALSE)
for (i in 1:length(spec)){
for (j in 1:length(measures)){
model_exp= sentiment_reg("dji_ret", measures[j], spec[i], "exp", "usrecd", sentiment_data)
model_rec= sentiment_reg("dji_ret", measures[j], spec[i], "rec", "usrecd", sentiment_data)
r2_exp= summary(model_exp)$r.squared
r2_rec= summary(model_rec)$r.squared
adjr2_exp= summary(model_exp)$adj.r.squared
adjr2_rec= summary(model_rec)$adj.r.squared
col_adjr2_exp= paste0("adjr2_exp_", spec[i])
col_adjr2_rec= paste0("adjr2_rec_", spec[i])
col_r2_exp= paste0("r2_exp_", spec[i])
col_r2_rec= paste0("r2_rec_", spec[i])
out_stats_r2[j, col_adjr2_exp]= adjr2_exp
out_stats_r2[j, col_adjr2_rec]= adjr2_rec
out_stats_r2[j, col_r2_exp]= r2_exp
out_stats_r2[j, col_r2_rec]= r2_rec
}
}
#Table 2
out_stats_r2_scaled= out_stats_r2
out_stats_r2_scaled[, -1]= out_stats_r2_scaled[, -1] * 100
stargazer(out_stats_r2_scaled, summary=F, digits= 2, rownames=F)
#Create lags
create_lags= function(data, column, n_lags) {
lag_list= setNames(
lapply(1:n_lags, function(i) {
dplyr::lag(data[[column]], i)
}),
paste0("lag", 1:n_lags, "_", column)
)
data %>%
mutate(!!!lag_list)
}
run_sentiment_simulations= function(nsim, data, dep_var, sentiment_var= "noise",
specification= "full", recession_ind= "usrecd",
n_lags= 5, noise_mean= 0, noise_sd= 1)
{
sim_results= vector("list", nsim)
for (k in seq_len(nsim)) {
data_sim= data
#Simulate noise
data_sim[[sentiment_var]]= rnorm(nrow(data_sim), mean= noise_mean, sd= noise_sd)
#Create lags
data_sim= create_lags(data_sim, sentiment_var, n_lags)
#Run predictive regressions
model_exp= sentiment_reg(dep_var, sentiment_var, specification, "exp",
recession_ind, data_sim)
model_rec= sentiment_reg(dep_var, sentiment_var, specification, "rec",
recession_ind, data_sim)
#Store betas and residuals
sim_results[[k]]= list(expansion = list(beta= coef(model_exp)[sentiment_var],
residuals= resid(model_exp)),
recession = list(beta= coef(model_rec)[sentiment_var],
residuals= resid(model_rec)))
}
return(sim_results)
}
#Uncomment lines below to run simulations and store them.
#Commented out for speed; running and storing takes ~10 minutes on a base M4 Mac Mini.
# sim_results_dji_ret= run_sentiment_simulations(nsim= 10000, data= sentiment_data,
#                                        dep_var= "dji_ret", sentiment_var= "noise",
#                                        specification= "full", recession_ind= "usrecd",
#                                        n_lags= 5)
#
# sim_results_dji_novol= run_sentiment_simulations(nsim= 10000, data= sentiment_data,
#                                                dep_var= "dji_novol", sentiment_var= "noise",
#                                                specification= "full", recession_ind= "usrecd",
#                                                n_lags= 5)
# saveRDS(sim_results_dji_ret, "sim_results_dji_ret.rds")
# saveRDS(sim_results_dji_novol, "sim_results_dji_novol.rds")
sim_results_dji_ret= readRDS("Simulations/sim_results_dji_ret.rds")
sim_results_dji_novol= readRDS("Simulations/sim_results_dji_novol.rds")
sim_lists= list(dji_ret= sim_results_dji_ret,  dji_novol= sim_results_dji_novol)
summary_table_sims= data.frame(series= character(),
mae_exp_mean= numeric(), mae_rec_mean= numeric(),
mse_exp_mean= numeric(), mse_rec_mean= numeric(),
beta_ratio_mean= numeric(), beta_ratio_mean_winsor= numeric(),
stringsAsFactors= FALSE)
for(s in names(sim_lists)){
sim_results= sim_lists[[s]]
nsim= length(sim_results)
mae_exp= numeric(nsim)
mae_rec= numeric(nsim)
mse_exp= numeric(nsim)
mse_rec= numeric(nsim)
beta_ratio= numeric(nsim)
for(k in seq_len(nsim)){
exp_res= sim_results[[k]]$expansion$residuals
rec_res= sim_results[[k]]$recession$residuals
mae_exp[k]= mean(abs(exp_res), na.rm= TRUE)
mae_rec[k]= mean(abs(rec_res), na.rm= TRUE)
mse_exp[k]= mean(exp_res^2, na.rm= TRUE)
mse_rec[k]= mean(rec_res^2, na.rm= TRUE)
beta_exp= sim_results[[k]]$expansion$beta
beta_rec= sim_results[[k]]$recession$beta
beta_ratio[k]= abs(beta_rec / beta_exp)
}
if(s == "dji_ret") beta_ratios_dji_ret= beta_ratio
if(s == "dji_novol") beta_ratios_dji_novol= beta_ratio
summary_table_sims= rbind(summary_table_sims,
data.frame(series= s,
mae_exp_mean= mean(mae_exp),
mse_exp_mean= mean(mse_exp),
mae_rec_mean= mean(mae_rec),
mse_rec_mean= mean(mse_rec),
beta_ratio_mean= mean(beta_ratio),
beta_ratio_mean_winsor= mean(winsorize(beta_ratio,
threshold= 0.01, method= "percentile")),
stringsAsFactors= FALSE))
}
winsor_beta_ratio_dji_ret= winsorize(beta_ratios_dji_ret, threshold= 0.01, method= "percentile")
winsor_beta_ratio_dji_novol= winsorize(beta_ratios_dji_novol, threshold= 0.01, method= "percentile")
make_binned_df= function(x, series_name){
bins= cut(x,
breaks= c(seq(0, 10, 1), Inf),
include.lowest= TRUE,
right= FALSE)
levels(bins)= c(
paste0("[", 0:9, "–", 1:10, ")"),
"[10–∞)"
)
data.frame(
value= x,
bin= bins,
series= series_name
)
}
plot_beta_hist= function(df, vlines){
df$bin_num= as.numeric(df$bin)
bin_labels= levels(df$bin)
mean_val= mean(df$value, na.rm= TRUE)
print(mean_val)
ggplot(df, aes(x= bin_num)) +
geom_bar(fill= "#006699", color= "white") +
# vertical dashed reference lines
geom_vline(
xintercept= vlines,
linetype= "dashed",
linewidth= 2.5,
color= "black"
) +
# mean line (solid)
geom_vline(
xintercept= mean_val+0.5,
linetype= "solid",
linewidth= 2.5,
color= "black"
) +
scale_x_continuous(
breaks= seq_along(bin_labels),
labels= bin_labels
) +
theme_minimal() +
labs(x= "Beta Ratio", y= "Frequency") +
theme(
axis.title= element_text(size= 72, face= "bold"),
axis.text= element_text(size= 72, face= "bold"),
axis.text.x= element_text(
size= 72,
face= "bold",
angle= 45,
hjust= 1,
vjust= 1))
}
#Figure 1, save with width 2400 and maintain aspect ratio
df_ret= make_binned_df(winsor_beta_ratio_dji_ret, "dji_ret")
df_novol= make_binned_df(winsor_beta_ratio_dji_novol, "dji_novol")
#Add 0.5 so the lines appear correct on the plot (columns are centered at .5)
vlines_ret= c(abs(out_stats_beta_scaled[1,3]/out_stats_beta_scaled[1,2])+0.5,
abs(out_stats_beta_scaled[1,7]/out_stats_beta_scaled[1,6])+0.5,
abs(out_stats_beta_scaled[1,11]/out_stats_beta_scaled[1,10])+0.5)
vlines_novol= c(abs(out_stats_beta_scaled[2,3]/out_stats_beta_scaled[2,2])+0.5,
abs(out_stats_beta_scaled[2,7]/out_stats_beta_scaled[2,6])+0.5,
abs(out_stats_beta_scaled[2,11]/out_stats_beta_scaled[2,10])+0.5)
plot_beta_hist(df_ret, vlines_ret)
plot_beta_hist(df_novol, vlines_novol)
#Table 7
summary_table_sims_scaled= summary_table_sims
summary_table_sims_scaled[1, c(2, 4)]= summary_table_sims_scaled[1, c(2, 4)] * 10000
summary_table_sims_scaled[1, c(3, 5)]= summary_table_sims_scaled[1, c(3, 5)] * 100000000
summary_table_sims_scaled[2, c(2, 4)]= summary_table_sims_scaled[2, c(2, 4)] * 100
summary_table_sims_scaled[2, c(3, 5)]= summary_table_sims_scaled[2, c(3, 5)] * 10000
stargazer(summary_table_sims_scaled, summary=F, digits= 2)
avail_series= c("dji_ret", "dji_novol")
measures= c("measure_pos_unit", "measure_neg_unit", "measure_pess_unit")
sim_lists= list(dji_ret= sim_results_dji_ret,
dji_novol= sim_results_dji_novol)
dm_results1= data.frame(measure= character(), ret_series= character(),
dm_mae_exp1= numeric(), dm_mse_exp1= numeric(),
dm_mae_rec1= numeric(), dm_mse_rec1= numeric(),
stringsAsFactors= FALSE)
dm_results5= data.frame(measure= character(), ret_series= character(),
dm_mae_exp5= numeric(), dm_mse_exp5= numeric(),
dm_mae_rec5= numeric(), dm_mse_rec5= numeric(),
stringsAsFactors= FALSE)
oos_r2_results= data.frame(measure= character(), ret_series= character(),
oos_r2_exp_mean= numeric(), oos_r2_exp_sd= numeric(),
oos_r2_exp_min= numeric(), oos_r2_exp_max= numeric(),
oos_r2_rec_mean= numeric(), oos_r2_rec_sd= numeric(),
oos_r2_rec_min= numeric(), oos_r2_rec_max= numeric(),
stringsAsFactors = FALSE)
calc_r2_oos = function(real, bench){
real = real[!is.na(real)]
bench = bench[!is.na(bench)]
1 - mean(real^2) / mean(bench^2)
}
for (i in seq_along(avail_series)) {
current_series= avail_series[i]
sim_results= sim_lists[[current_series]]
for (j in seq_along(measures)) {
current_measure= measures[j]
model_exp= sentiment_reg(current_series, current_measure, "full",
"exp", "usrecd", sentiment_data)
model_rec= sentiment_reg(current_series, current_measure, "full",
"rec", "usrecd", sentiment_data)
actual_resid_exp= resid(model_exp)
actual_resid_rec= resid(model_rec)
nsim= length(sim_results)
oos_exp_vec= numeric(nsim)
oos_rec_vec= numeric(nsim)
for (k in seq_len(nsim)) {
sim_resid_exp= sim_results[[k]]$expansion$residuals
sim_resid_rec= sim_results[[k]]$recession$residuals
oos_exp_vec[k]= calc_r2_oos(actual_resid_exp, sim_resid_exp)
oos_rec_vec[k]= calc_r2_oos(actual_resid_rec, sim_resid_rec)
}
clean_measure= gsub("measure_|_unit", "", current_measure)
oos_r2_results= rbind(oos_r2_results,
data.frame(measure= clean_measure,
ret_series= current_series,
oos_r2_exp_mean= mean(oos_exp_vec, na.rm= TRUE),
oos_r2_exp_sd = sd(oos_exp_vec, na.rm= TRUE),
oos_r2_exp_min = min(oos_exp_vec, na.rm= TRUE),
oos_r2_exp_max = max(oos_exp_vec, na.rm= TRUE),
oos_r2_rec_mean = mean(oos_rec_vec, na.rm= TRUE),
oos_r2_rec_sd = sd(oos_rec_vec, na.rm= TRUE),
oos_r2_rec_min = min(oos_rec_vec, na.rm= TRUE),
oos_r2_rec_max = max(oos_rec_vec, na.rm= TRUE),
stringsAsFactors = FALSE))
}
}
oos_tab= oos_r2_results
oos_tab[ , sapply(oos_tab, is.numeric)]= oos_tab[ , sapply(oos_tab, is.numeric)] * 100
stargazer(oos_tab, summary=F, digits= 3, rownames=F)
View(oos_r2_results)
